{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rzn3xV2ScXXN",
        "papermill": {
          "duration": 0.013538,
          "end_time": "2021-02-20T13:42:41.189745",
          "exception": false,
          "start_time": "2021-02-20T13:42:41.176207",
          "status": "completed"
        },
        "tags": []
      },
      "source": [
        "# Analyzing RCT data with Precision Adjustment"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "StJRPl8Lcsdd"
      },
      "outputs": [],
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "import patsy\n",
        "import statsmodels.formula.api as smf\n",
        "import statsmodels.api as sm"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "U1aS4EabcXXr",
        "papermill": {
          "duration": 0.011347,
          "end_time": "2021-02-20T13:42:41.213041",
          "exception": false,
          "start_time": "2021-02-20T13:42:41.201694",
          "status": "completed"
        },
        "tags": []
      },
      "source": [
        "## Data\n",
        "\n",
        "In this lab, we analyze the Pennsylvania re-employment bonus experiment, which was previously studied in \"Sequential testing of duration data: the case of the Pennsylvania ‘reemployment bonus’ experiment\" (Bilias, 2000), among others. These experiments were conducted in the 1980s by the U.S. Department of Labor to test the incentive effects of alternative compensation schemes for unemployment insurance (UI). In these experiments, UI claimants were randomly assigned either to a control group or one of five treatment groups. Actually, there are six treatment groups in the experiments. Here we focus on treatment group 4, but feel free to explore other treatment groups. In the control group the current rules of the UI applied. Individuals in the treatment groups were offered a cash bonus if they found a job within some pre-specified period of time (qualification period), provided that the job was retained for a specified duration. The treatments differed in the level of the bonus, the length of the qualification period, and whether the bonus was declining over time in the qualification period; see http://qed.econ.queensu.ca/jae/2000-v15.6/bilias/readme.b.txt for further details on data. \n",
        "  "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "2qKDD0hecXX4",
        "papermill": {
          "duration": 0.279482,
          "end_time": "2021-02-20T13:42:41.503829",
          "exception": false,
          "start_time": "2021-02-20T13:42:41.224347",
          "status": "completed"
        },
        "tags": []
      },
      "outputs": [],
      "source": [
        "data = pd.read_csv(\"https://raw.githubusercontent.com/VC2015/DMLonGitHub/master/penn_jae.dat\", delim_whitespace=True)\n",
        "n, p = data.shape\n",
        "data = data[data[\"tg\"].isin([0, 4])]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CWoXcstjcXX6",
        "papermill": {
          "duration": 0.065266,
          "end_time": "2021-02-20T13:42:41.580972",
          "exception": false,
          "start_time": "2021-02-20T13:42:41.515706",
          "status": "completed"
        },
        "tags": []
      },
      "outputs": [],
      "source": [
        "data[\"T4\"] = np.where(data[\"tg\"] == 4, 1, 0)\n",
        "data[\"T4\"].value_counts()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5q2hZMxScXX8",
        "papermill": {
          "duration": 0.053285,
          "end_time": "2021-02-20T13:42:41.646704",
          "exception": false,
          "start_time": "2021-02-20T13:42:41.593419",
          "status": "completed"
        },
        "tags": []
      },
      "outputs": [],
      "source": [
        "data.describe()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bgCUMjFWcXYB",
        "papermill": {
          "duration": 0.012846,
          "end_time": "2021-02-20T13:42:41.672662",
          "exception": false,
          "start_time": "2021-02-20T13:42:41.659816",
          "status": "completed"
        },
        "tags": []
      },
      "source": [
        "### Model \n",
        "To evaluate the impact of the treatments on unemployment duration, we consider the linear regression model:\n",
        "\n",
        "$$\n",
        "Y =  D \\beta_1 + W'\\beta_2 + \\varepsilon, \\quad E \\varepsilon (D,W')' = 0,\n",
        "$$\n",
        "\n",
        "where $Y$ is  the  log of duration of unemployment, $D$ is a treatment  indicators,  and $W$ is a set of controls including age group dummies, gender, race, number of dependents, quarter of the experiment, location within the state, existence of recall expectations, and type of occupation.   Here $\\beta_1$ is the ATE, if the RCT assumptions hold rigorously.\n",
        "\n",
        "\n",
        "We also consider interactive regression model:\n",
        "\n",
        "$$\n",
        "Y =  D \\alpha_1 + D W' \\alpha_2 + W'\\beta_2 + \\varepsilon, \\quad E \\varepsilon (D,W', DW')' = 0,\n",
        "$$\n",
        "where $W$'s are demeaned (apart from the intercept), so that $\\alpha_1$ is the ATE, if the RCT assumptions hold rigorously."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CfgwXd6ycXYC",
        "papermill": {
          "duration": 0.012735,
          "end_time": "2021-02-20T13:42:41.698676",
          "exception": false,
          "start_time": "2021-02-20T13:42:41.685941",
          "status": "completed"
        },
        "tags": []
      },
      "source": [
        "Under RCT, the projection coefficient $\\beta_1$ has\n",
        "the interpretation of the causal effect of the treatment on\n",
        "the average outcome. We thus refer to $\\beta_1$ as the average\n",
        "treatment effect (ATE). Note that the covariates, here are\n",
        "independent of the treatment $D$, so we can identify $\\beta_1$ by\n",
        "just linear regression of $Y$ on $D$, without adding covariates.\n",
        "However we do add covariates in an effort to improve the\n",
        "precision of our estimates of the average treatment effect."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XXlBcN0ucXYC",
        "papermill": {
          "duration": 0.012728,
          "end_time": "2021-02-20T13:42:41.724182",
          "exception": false,
          "start_time": "2021-02-20T13:42:41.711454",
          "status": "completed"
        },
        "tags": []
      },
      "source": [
        "### Analysis\n",
        "\n",
        "We consider \n",
        "\n",
        "*  classical 2-sample approach, no adjustment (CL)\n",
        "*  classical linear regression adjustment (CRA)\n",
        "*  interactive regression adjusment (IRA)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pWqd1FIJcXYD",
        "papermill": {
          "duration": 0.01276,
          "end_time": "2021-02-20T13:42:41.749736",
          "exception": false,
          "start_time": "2021-02-20T13:42:41.736976",
          "status": "completed"
        },
        "tags": []
      },
      "source": [
        "# Carry out covariate balance check"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "A3GrR6ZlcXYE",
        "papermill": {
          "duration": 0.013486,
          "end_time": "2021-02-20T13:42:41.776684",
          "exception": false,
          "start_time": "2021-02-20T13:42:41.763198",
          "status": "completed"
        },
        "tags": []
      },
      "source": [
        "This is done using \"cov_type='HC1'\" which computes the correct Eicher-Huber-White standard errors, instead of the classical non-robust formula based on the homoscdedasticity."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Jnh2fbqRcXYE",
        "papermill": {
          "duration": 0.428084,
          "end_time": "2021-02-20T13:42:42.217973",
          "exception": false,
          "start_time": "2021-02-20T13:42:41.789889",
          "status": "completed"
        },
        "tags": []
      },
      "outputs": [],
      "source": [
        "covariate_balance = smf.ols(\"T4~(female+black+othrace+C(dep)+q2+q3+q4+q5+q6+agelt35+agegt54+durable+lusd+husd)**2\",\n",
        "                            data)\n",
        "covariate_balance.fit(cov_type=\"HC1\").summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0AC_UNXGcXYF",
        "papermill": {
          "duration": 0.02549,
          "end_time": "2021-02-20T13:42:42.269256",
          "exception": false,
          "start_time": "2021-02-20T13:42:42.243766",
          "status": "completed"
        },
        "tags": []
      },
      "source": [
        "We see that that even though this is a randomized experiment, balance conditions are failed."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AlZe_Dr6cXYH",
        "papermill": {
          "duration": 0.021046,
          "end_time": "2021-02-20T13:42:42.315674",
          "exception": false,
          "start_time": "2021-02-20T13:42:42.294628",
          "status": "completed"
        },
        "tags": []
      },
      "source": [
        "# Model Specification"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CQPB5YfCcXYH",
        "papermill": {
          "duration": 0.350081,
          "end_time": "2021-02-20T13:42:42.680497",
          "exception": false,
          "start_time": "2021-02-20T13:42:42.330416",
          "status": "completed"
        },
        "tags": []
      },
      "outputs": [],
      "source": [
        "# no adjustment (2-sample approach)\n",
        "cl = smf.ols(\"np.log(inuidur1) ~ T4\", data)\n",
        "\n",
        "# adding controls\n",
        "cra = smf.ols(\"np.log(inuidur1) ~ T4 + (female+black+othrace+C(dep)+q2+q3+q4+q5+q6+agelt35+agegt54+durable+lusd+husd)**2\",\n",
        "              data)\n",
        "# Omitted dummies: q1, nondurable, muld\n",
        "\n",
        "cl_results = cl.fit(cov_type=\"HC1\")\n",
        "cl_results.summary()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "EQl_YGSLNWAi"
      },
      "outputs": [],
      "source": [
        "cra_results = cra.fit(cov_type=\"HC1\")\n",
        "cra_results.summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "R-G1Vqo7cXYI",
        "papermill": {
          "duration": 0.02768,
          "end_time": "2021-02-20T13:42:42.735966",
          "exception": false,
          "start_time": "2021-02-20T13:42:42.708286",
          "status": "completed"
        },
        "tags": []
      },
      "source": [
        "The interactive specification corresponds to the approach introduced in Lin (2013)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "SWBol4-GcXYJ",
        "papermill": {
          "duration": 0.898452,
          "end_time": "2021-02-20T13:42:43.662684",
          "exception": false,
          "start_time": "2021-02-20T13:42:42.764232",
          "status": "completed"
        },
        "tags": []
      },
      "outputs": [],
      "source": [
        "#interactive regression model\n",
        "\n",
        "ira_formula = \"0 + (female+black+othrace+C(dep)+q2+q3+q4+q5+q6+agelt35+agegt54+durable+lusd+husd)**2\"\n",
        "X = patsy.dmatrix(ira_formula, data, return_type='dataframe')\n",
        "X.columns = [f'x{t}' for t in range(X.shape[1])] # clean column names\n",
        "X = (X - X.mean(axis=0)) # demean all control covariates\n",
        "\n",
        "# construct interactions of treatment and (de-meaned covariates, 1)\n",
        "ira_formula = \"T4 * (\"+ \"+\".join(X.columns) + \")\"\n",
        "X['T4'] = data['T4']\n",
        "X = patsy.dmatrix(ira_formula, X, return_type='dataframe')\n",
        "ira = sm.OLS(np.log(data[[\"inuidur1\"]]), X)\n",
        "ira_results = ira.fit(cov_type=\"HC1\")\n",
        "ira_results.summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ovx__MnpcXYK",
        "papermill": {
          "duration": 0.030243,
          "end_time": "2021-02-20T13:42:43.724231",
          "exception": false,
          "start_time": "2021-02-20T13:42:43.693988",
          "status": "completed"
        },
        "tags": []
      },
      "source": [
        "Next we try out partialling out with lasso"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "eWcw5TBncXYK",
        "papermill": {
          "duration": 1.29547,
          "end_time": "2021-02-20T13:42:45.045172",
          "exception": false,
          "start_time": "2021-02-20T13:42:43.749702",
          "status": "completed"
        },
        "tags": []
      },
      "outputs": [],
      "source": [
        "from sklearn.linear_model import LassoCV\n",
        "\n",
        "D = X['T4']\n",
        "A = X.drop('T4', axis=1) # contains X and interaction terms D*X\n",
        "y = np.log(data[\"inuidur1\"])\n",
        "\n",
        "Dres = D - np.mean(D)\n",
        "yres = y - LassoCV(max_iter=10000).fit(A, y).predict(A)\n",
        "\n",
        "ira_lasso_results = sm.OLS(yres, Dres).fit(cov_type='HC1')\n",
        "ira_lasso_results.summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TNAQW7KwcXYQ",
        "papermill": {
          "duration": 0.030965,
          "end_time": "2021-02-20T13:42:45.107753",
          "exception": false,
          "start_time": "2021-02-20T13:42:45.076788",
          "status": "completed"
        },
        "tags": []
      },
      "source": [
        "### Results"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "RmPEWZXxcXYT",
        "papermill": {
          "duration": 0.113372,
          "end_time": "2021-02-20T13:42:45.304786",
          "exception": false,
          "start_time": "2021-02-20T13:42:45.191414",
          "status": "completed"
        },
        "tags": []
      },
      "outputs": [],
      "source": [
        "results = {\n",
        "    \"CL\": [cl_results.params[cl.exog_names.index(\"T4\")], \n",
        "           cl_results.HC1_se[cl.exog_names.index(\"T4\")]],\n",
        "    \"CRA\": [cra_results.params[cra.exog_names.index(\"T4\")], \n",
        "           cra_results.HC1_se[cra.exog_names.index(\"T4\")]],\n",
        "    \"IRA\": [ira_results.params[ira.exog_names.index(\"T4\")], \n",
        "           ira_results.HC1_se[ira.exog_names.index(\"T4\")]],\n",
        "    \"IRA w Lasso\": [ira_lasso_results.params[0], \n",
        "           ira_lasso_results.bse[0]],\n",
        "}\n",
        "results = pd.DataFrame(results)\n",
        "results.index = [\"Estimate\", \"Standard Error\"]\n",
        "results\n",
        "# for printing to LaTeX: print(results.to_latex())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tT1XCKLhcXYU",
        "papermill": {
          "duration": 0.019862,
          "end_time": "2021-02-20T13:42:45.344887",
          "exception": false,
          "start_time": "2021-02-20T13:42:45.325025",
          "status": "completed"
        },
        "tags": []
      },
      "source": [
        "Treatment group 4 experiences an average decrease of about $7.8\\%$ in the length of unemployment spell.\n",
        "\n",
        "\n",
        "Observe that regression estimators delivers estimates that are slighly more efficient (lower standard errors) than the simple 2 mean estimator, but essentially all methods have very similar standard errors. From IRA results we also see that there is not any statistically detectable heterogeneity.  We also see the regression estimators offer slightly lower estimates -- these difference occur perhaps to due minor imbalance in the treatment allocation, which the regression estimators try to correct.\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "odb5eio_AUfI"
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "colab": {
      "name": "analyzing-rct-reemployment-experiment.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.13"
    },
    "papermill": {
      "default_parameters": {},
      "duration": 7.178818,
      "end_time": "2021-02-20T13:42:45.473549",
      "environment_variables": {},
      "exception": null,
      "input_path": "__notebook__.ipynb",
      "output_path": "__notebook__.ipynb",
      "parameters": {},
      "start_time": "2021-02-20T13:42:38.294731",
      "version": "2.2.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}
